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Abstract 

This paper deals with bounding the error on the estimation of quantities 
of interest obtained by finite element and domain decomposition methods. 

The proposed bounds are written in order to separate the two errors involved 
in the resolution of reference and adjoint problems : on the one hand the 
discretization error due to the finite element method and on the other hand 
the algebraic error due to the use of the iterative solver. Beside practical 
considerations on the parallel computation of the bounds, it is shown that 
the interface conformity can be slightly relaxed so that local enrichment or 
refinement are possible in the subdomains bearing singularities or quantities 
of interest which simplifies the improvement of the estimation. Academic 
assessments are given on 2D static linear mechanic problems. 

^NOTICE: this is the author’s version of a work that was accepted for publication in CMAME. 
Changes resulting from the publishing process, such as peer review, editing, corrections, struc¬ 
tural formatting, and other quality control mechanisms may not be reflected in this document. 
Changes may have been made to this work since it was submitted for publication. A definitive 
version was subsequently published in Computer methods in Applied Mechanics and Engineering, 
2015, doi: 10.1016/j.cma.2015.01.009. 
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1 Introduction 

In order to certify structures by computations, two fundamental abilities must 
be gathered: (i) conducting large computations which can be achieved by using 
non-overlapping domain decomposition methods to exploit parallel hardware PZH 
mm, (n) warrantying the quality of the results which implies to use verification 
techniques. The computed approximation is then completed by an estimation of 
its distance to the exact solution, that distance being either on a global norm 
m na uni ED] or on a selection of engineering quantities of interest [311122112]. 

This paper is the sequel of contributions by the authors in order to associate 
domain decomposition methods (DDM) and verification techniques. In [28] it was 
shown how admissible fields (aka balanced residuals) could be naturally computed 
in parallel when solving one finite element system with classical DDM solvers 
like FETI or BDD. This tool was necessary to obtain guaranteed error bounds. 
Moreover the proposed method encapsulated classical sequential techniques [181 
EolEZlEg that could be employed independently on the subdomains. In EH a 
global upper bound of the error was separated in two: the first contribution was 
purely due to the iterative solver employed in DDM and it could be made as small 
as wanted by doing more iterations, whereas the second contribution depended 
on the discretization and it was quasi-const ant during the iterations. This result 
enabled us to stop iterations when the overall quality of the approximation could 
not be improved, which was much earlier than classical criteria on the norm of the 
residual would have suggested. 

The purpose of this paper is to extend these works to the handling of quantities 
of interest. First we show how estimation on quantities of interest can be conducted 
easily within DDM with a separation of the contribution of the solver and of 
the discretization. Second we demonstrate that all the wanted properties can be 
preserved with nested kinematics at the interface which simplifies the improvement 
of the estimation by local adaptation. 

The question of separating the contribution of the solver and of the discretiza¬ 
tion, in goal-oriented error estimation, was addressed in several other studies with 
different approaches. In [251124] , the quantity of interest was the mean value of the 
unknown held in a region of the structure, the coarse mesh of the structure was 
seen as a decomposition of the true refined mesh which allowed to take advantage 
of the properties of FETI algorithm. In [2T5], the authors proposed constant-free 
asymptotic upper and lower bounds but they did not connect the iteration and 
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discretization contributions. In HU. error indicators separating contributions were 
developed in the context of multigrids methods for Poisson and Stockes equations. 
They lead to a gain in terms of CPU time but the bounds were not guaranteed. 
Using the dual weighted residual method [2j for goal-oriented error estimation, a 
balance of discretization and iteration errors was proposed for eigenvalue prob¬ 
lems in [32] . Finally, in [40], the authors proposed goal-oriented a posteriori error 
estimator for problems discretized with mixed finite elements and multiscale do¬ 
main decomposition that separated the contributions of the discretization of the 
subdomains from the contribution of the mortar interface. 

In this paper, we employ classical extractor techniques [251 39j 22], which leads 
to the definition of an adjoint problem which is solved at the same time as the for¬ 
ward problem using a block Krylov algorithm |£|. Thanks to [281 [34] we compute 
guaranteed upper and lower bounds of the quantity of interest at each iteration 
with the separation of the two sources of error for each problem (forward and ad¬ 
joint), which enables us to drive the iterative solver by the error on the quantity 
of interest. 

Moreover, since we consider quantities of interest defined on small supports, 
the adjoint problems have very localized solution. A simple strategy to improve 
the quality of the estimation is then to better the resolution of the adjoint prob¬ 
lem either with a locally refined mesh [5] or using the partition of unity method 
and handbook techniques pjj. We show that within the framework of domain de¬ 
composition methods, it is possible to enrich or refine the kinematic of selected 
subdomains without impairing the exactness of the bound nor the ease of compu¬ 
tation even if a certain incompatibility appears at the interface. 

The paper is organized as follows. In Section [2] we present the goal-oriented 
error estimation within domain decomposition framework; we propose bounds on 
quantities of interest that separate the contributions of the discretization and of 
the solver. In Section [3] we show that nested interfaces can be employed which 
simplifies the improvement of the quality of the bounds by local refinement. As¬ 
sessments are presented in Section [4] and Section [5] concludes the paper. 
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2 Error estimation on quantity of interest in sub¬ 
structured context 

2.1 Substructured formulation of the forward and adjoint 
problems 

2.1.1 Reference problems 

Let represent the physical space (d = 2 or 3). Let us consider the static 
equilibrium of a (polyhedral) structure which occupies the open domain 12 c W l 
and which is subjected to given body force / within 12, to given traction force 
g on d g Q and to given displacement held u d on the complementary part of the 
boundary d u Q (meas(d u 12) ^ 0). We assume that the structure undergoes small 
perturbations and that the material is linear elastic, characterized by Hooke’s 
elasticity tensor EL Let u be the unknown displacement held, e (u) the symmetric 
part of the gradient of u, a the Cauchy stress tensor. In order to evaluate (localized) 
quantities of interest, we consider an adjoint problem loaded by a (linear) extractor 
l with small support. Fields related to the adjoint problem will always be noted 
with a tilde (for instance adjoint displacement is u). 

Let u c 12 be an open subset of 12. We introduce three affine subspaces, one 
linear space and one positive form: 

• Affine subspace of kinematic admissible helds (KA-helds) 

KA(w) = ju e (H 1 (o;)) d , u = Ud on dui d u 12 j (1) 

• Associated linear subspace KA°: 

KA(w) = KA°(u;) = jn e (H 1 (o;)) <i , u = 0 on da;\d 9 Q| (2) 


Affine subspace of statically admissible helds (SA-helds) for the forward prob¬ 
lem: 


SA(w) = \ L e (L 2 M)^; m d ; Vn e KA°( W ), 


r:e(v)dQ = I / • vdQ + J g ■ vdS = l(v) )• (3) 

d g O P| dui 


We note l the linear form which gathers the loads of the forward problem. 
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• Affine subspace of statically admissible fields (SA-fields) for the adjoint prob¬ 
lem: 

SA(w) = jr e (L 2 (w))^; Vn e KA 0 (w), jL ■ g (v) du; = l(v) 1 (4) 

^ UJ ' 

• Error in constitutive equation 

ecRu (u, q) = \\q - H : g (u) || H -v (5) 

where 

The forward and adjoint problem set on fl can be formulated as: 

Find (v^^qjj e KA(Q) x SA(Q) such that e CRn {u ex ,q ex ) = 0 (6) 

Find e KA°(0) x SA(fi) such that e CRn (u ex ,q ex ) = 0 (7) 

The solution to these problems, named “exact” solutions, exist and are unique. 



2.1.2 Substructured formulation 

Let us consider a decomposition of domain SI in N s d open subsets such that 
P| ) = 0 for s A s' and Q = P s . The interface between subdomains 
r LA) = q ( s ) p| q(«') 

is supposed to be regular enough for traces of locally admissible 
fields to be well defined, n is the outer normal vector. 

The mechanical problem on the substructured configuration writes : 


' u {s) e KA(fi (s) ) 

Vs < q^ e SA(Q^ s) ) and V(s, s') 

. ecR nW (uW,gW) = 0 


tr(u (s) ) = tr(u (s,) ) on T (s ’ s,) 

a {s) . n (s) + a {s') . n {s') = o on r (s,0 ^ 


We assume the same domain decomposition is used for the forward and adjoint 
problems. The adjoint problem satisfies the same type of substructured formu¬ 
lation: the fields u (respectively u) and q (resp. q) are globally admissible if 
(w (s) ,a (s) ) (resp. (u^ s \ q ^)) are kinematically and statically admissible inside each 
and if the displacements are continuous and tractions balanced on interfaces. 
The set of fields u (or u) dehned on SI such that they are kinematically admissible 
inside each domain without interface continuity is a broken space which we 
note KA(UfiW) and KA°((Jf2^). Similarly, we denote SA(|JSA S )) the broken 
space of statically admissible helds inside subdomains. 
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Remark. The error in constitutive relation is by essence a parallel quantity: 

V«6KA(Un w ),V£eSA(Un<*>), e CHa (u,g) = ^2>cs n( . l (<‘ ( * ) .g W ) (9) 

2.2 Finite element approximation for the forward and ad¬ 
joint problems 

Both forward and adjoint problems are solved with a finite element approxima¬ 
tion. Note that the mesh used to solve the adjoint problem can be completely 
different from the mesh used for the forward problem. We detail the finite element 
approximation for the forward problem (the principle is the same for the adjoint 
problem). 

2.2.1 Global form 

We associate with the mesh of fi the finite-element subspace KA#(12). The finite 
element approximation consists in searching: 

u H e KA h (Q), q_ R = H : £ (u H ) 

f g H '■ g(v H ) dtt = f / • v H (in + f g ■ v H dS = l(v H ), dv H e KA^(fi) 

Jf2 JO Jdnfl 

(io) 

We introduce the matrix ip H of shape functions which form a basis of KA#(f2) 
and the vector of nodal unknowns u. Therefore, we have u H = and obtain 
the following linear system: 

Ku = f (11) 

where K is the (symmetric positive definite) stiffness matrix and f is the vector of 
generalized forces (which includes imposed displacement). 

2.2.2 Substructured form 

For both problems, we assume that the mesh of 11 and the substructuring are con¬ 
forming. This hypothesis implies that each element only belongs to one subdomain 
and nodes are matching on the interfaces. Each degree of freedom is either located 
inside a subdomain (Subscript i) or on its boundary (Subscript b). 

Let be the discrete trace operator, so that u^ = t^u^). Let us intro¬ 
duce the unknown nodal reaction on the interface the equilibrium of each 
subdomain writes: 

K (s) u (s) = f (s) + t (s)T A (s) (12) 
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Let (A*- 5 -*) and (B 1 ^) be the primal and dual assembly operator so that the discrete 
counterpart of the interface admissibility equations is: 


VB (s) t (s) u (s) = 0 


S 


(13) 



A (,) A<*> = 0 


Note that in the case when Subdomain s has not enough Dirichlet boundary 


conditions then the reaction has to balance the subdomain with respect to rigid 
body motions. Let be a basis of ker(K^ s ^), we have: 


R( s)T (f (s ) + t (s)T A (s) ) = 0 


(14) 


2.2.3 Domain decompostion solvers 

First, let us briefly present the two classical domain decomposition solvers BDD 
and FETI (see Algorithm [2] with for now M = N = I for BDD and [34] for 
FETI, see also Pi and its bibliography for more details). In BDD solver, a 
unique interface displacement is introduced so that the continuity of displacement 
is automatically verified. Using the primal Schur complement, the problem can be 
rewritten with interface quantities. In FETI solver, a unique set of nodal forces 
on the interface is introduced so that the balance of forces is verified. The dual 
Schur complement allows to write the problem with interface quantities. 

Then, the resolution of interface problem is based on iterative Krylov solver, 
namely a preconditioned conjugate gradient. The classical preconditioners make 
use of scaled assembling operators (A*®)) and (Eh s )) solutions of equations [35, 12] : 


Yj A 1 *'A ( ' )T - I and V B w B ( * )r - I 


(15) 


s s 


During the algorithm, two types of mechanical problems are solved locally (on 


each subdomain): problems with Dirichlet conditions on the interface (subscript 
D) and problems with Neumann conditions on the interface (subscript N). They 
correspond to the following procedures: 
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When developing these methods, we get: 


(d?)i = (f<*> - Kg’u'*') and 

\( s ) _ c( s )ii( s ) f( s ) , tv( s )tv( s ) 1 f( s ) 

a. d — O u b I 6 -I- *^bi h 

u w = K (s)+ ( f (*) + t W T A^j 


( u d)& = u 


00 

b 


where 


K 


(s) 

bb 


*^bi rv ii **-ib 


is the well-known Schur complement ma¬ 


trix, and K^ + is a pseudo-inverse of K^. Thanks to projector (Px) and initial¬ 
ization (method Initialize), the Neumann problems are well-posed (see |10]). 


If the forward and adjoint problems are approximated on the same mesh, they 
share the same stiffness matrix and they can be solved simultaneously using a 
block-conjugate gradient. In Algorithm [2j we use capital letter to note blocks of 
forward/adjoint vectors: typically U = [u, u], same notation holds for the residual 
R = [r,r], preconditioned residual Z = [z,z], search directions W and Q; a is a 
2x2 matrix. One has to pay attention to the potential (yet unlikely) singularity 
of the 2x2 matrix (Q T W) which needs to be inverted. 

2.2.4 Description of admissible fields 

At every iteration of the solver, a totally parallel computation of the following 
fields is possible |28| : 

• u n = (y$) s e KA(fi) and u D = (v^) s e KA°(0): globally admissible 
displacement fields, (A^, A^) are the associated nodal reaction which are 
not balanced before convergence. 

• ILn = (Atv)s e KA(|Jf2^) and u N = (u$) s E KA°((JQ^): displacements 
associated to a given balanced Neumann condition A^ (resp. A^), they are 
not necessarily continuous across interfaces. 

• off = H : £ : th e s f ress field associated to u$ ■ It can be employed 

(with additional input A^ ) to build stress fields a_^ which are statically 
admissible g_ N = (<r^) s £ SA(fl). Same property holds for the adjoint 
problem, we can build o_n = (g^) s E SA(fi). 

As a consequence, for both forward and adjoint problems, kinematically displace¬ 
ments fields and statically stress fields are available even if the solver is not con¬ 
verged. Those fields are used to compute global error estimation for both problems. 


Indeed, we have the following guaranteed upper bounds which can be computed 
in parallel [28]: 


Ue X - Main ^ e CR n(u D ,g N ) 

Uex ^ eC Rn (u D , Ojy) 


(16) 


where j||x|||n = ||e(x) \\w~\n- 


2.3 Error estimation of quantities of interest 

As said earlier the adjoint problem is meant to extract one quantity of interest in 
the forward problem. Let I ex = l(u ex ) be the unknown exact value of the quantity 
of interest. l(u D ) is an approximation of this quantity of interest. 


2.3.1 Upper and lower bounds with separated contributions 

We demonstrate two results which provide upper and lower bounds on the quantity 
of interest, as a product of error estimates on the forward and adjoint problems 
with, for each of these factors, a separation of contributions of the solver and of 
the discretization. 


Property 1. Using notation of Algorithm^, we have 

\I e x - K^d) ~ fjfjf.il < (VUz + e CR n(u N ,g N ) S ) + e C R 0 (u N ,gN) > ) (17) 

with 

Ihh ,i = f (g N - EI| (u D )) : | (Md) dtt ( 18 ) 

Jn 

Proof. This result is a direct application of the results demonstrated in HH1I5]: 


| 4 z - KUld) ~ TtfJf.il ^ e G R n (u D ,^ Ar )eci? n (u D ,^v) 


(19) 


with Ihh, i given in (1 18 [) . The separation between the sources of error is obtained 
using the triangular inequality and Lemma 2 presented in [23] (which proves the 
equality \\\u N — u D |||^ = Vr T z), it slightly improves the result given in remark 25 
of |33]: 


ecR n {uD>g N ) = \\g N - H : |(wjj) lln-yn 

= Wgjf - El : I (ujv) + H | (u N - u D ) || H -i,n 
A — El : g (u N ) He-i.n + |||wjv — Holllsi 
< ecR n (u N ,g N ) + 


( 20 ) 


□ 
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The quantity Vr 7 z (respectively Vr T z) is a measure of the residual of the 
forward (resp. adjoint) problem: it corresponds to a pure algebraic error. The 
quantity ecR (a) (u N ,g_ N ) (resp. ecR n(s) (u N ,ojy))) is a parallel estimate of the dis¬ 
cretization error of the forward (resp. adjoint) problem. As a consequence, we 
can define a new stopping criterion for the solver: the idea is to stop iterations 
when the algebraic error becomes negligible in comparison with the discretization 
error. This new stopping criterion is described in alg. [I] 


Algorithm 1: DD block-solver with adapted stopping criterion 

Set a > 1, /3 5= 1, 5 > 1 and f3 ^ 1; 

Initialize, get (U^,A^,R, Z); 

Estimate discretization error for reference problem e 2 = e CR. 


and adjoint problem e 2 = ^ 


e 2 

s e CR. 


n( s ) 


(U 


A) -A) 
■n 


nM 


(“» 4w) 




while Vr 7 z > e/a and \ / r T z > e/a do 

while Vr T z > e/(a(3) and \/r r z > e/a/3 do 


Make Alg. [2] iterations, get (U^, A^ J ' 


end 

Update error estimators e 2 = Z s e C R. 
e 2 




2 (~(s) -Ah 

e CR n(s) > iyv y 


rA) 




and 


end 


The objective being to have the norm of the residual a times smaller than the 
estimation of the discretization error for both problems (a = S = 10 is a typical 
value), the coefficient (3 = (3 (typically 2) takes into account the small variation of 
the estimate e and e of the discretization error along iterations. 

A second (twice better) bound can be derived in a similar way using another 
result of mm 

Property 2. 

14b - 1(Md) - i hh,2 | < ^(Vr^z + e CRn (u N ,g N ))(v + e CRn (u N ,^f)) (21) 
with 

Ihh.2 = \ J (m d )J : HT 1 : + % (u D )^j dVt (22) 
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2.3.2 Properties of Ihh ,i 

If the forward and adjoint problems share the same mesh, the quantity Ihh ,i can 
be rewritten in terms of quantity on the interface, which simplifies its evaluation: 


Property 3. 

Ihh ,i = (23) 

S 

Proof. The forward problems verified by u D e KA(f2) and u ex e KA(O) with a 
testing held u D e KA°(Q) read: 




i (s) 


| (m?) : H : | («£’) an = 2 '“( 2 d) + £ A<f 

s s 

2 f g(«“): h : g (si?) an = 2 + £ f • «">) • a 

C J^S c c Jr« 


Since = (u$) s e KA U (12) and is the exact solution, we have: 

E J r(s) (A ex -A (s) )-^ ) ^ = 0 


while since <f N is statically admissible and («^) s e KA°(0), we have: 

Ihh, i = f (g^ - % (md)) : | (S D ) dQ = f - Mg (u D )) : g (S D ) dfi 

Jo Jo 

Therefore, by substracting equation (1241) from equation (1251) . we obtain (l23j) . 


(24) 

(25) 


(26) 


(27) 

□ 


When the solver has reached convergence, if the meshes of the forward and 
adjoint problems are identical, the term Ihha is equal to zero (as a consequence of 
the Galerkin orthogonality). Thanks to the new expression of Ihh, it we can prove 
the following property. 


Property 4. Using a BDD solver, at each iteration, Ihha = 0. 

Proof. Using the notation of BDD’s algorithm (alg. [2] with M = N = I), we have: 

Ihh, i = r T u (28) 


The iterations of the block-CG can be defined by the following Galerkin scheme m 


Uj e Uo + /Q 

R« -L 


( 29 ) 
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where U = [u, u] , R = [r, r] , and /Q is the Krylov subspace associated with the 
preconditioned operator and the initial block of residuals. This property implies 
that 

Ihh, i = r T u = Tq u 0 (30) 

In order to prove the nullity of this term, one has to look more in detail at BDD’s 
initialization. We refer to ra and simply recall that r 0 e range (P^) while u 0 e 
ker(Pi). □ 


2.4 Assessment 

Let us consider a rectangular linear elastic structure hi = [—3/; 3/] x [—3/; 3/] with 
homogeneous Dirichlet boundary conditions. Young modulus is 1 Pa and Poisson 
coefficient is 0.3. The domain is subjected to a polynomial body force such that 
the exact solution is known: 

Vex = ( x + 31)(x - 31)(y + 31)(y - 31) ((y - 3l) 2 e^ + (y + 3l)e y ) 

The quantity of interest is the average of the normal stress in the direction e y in 
a small square of side Z/10 centered at the origin. 

Ku) = -Vk f fey ® £y) : S (ft) du (31) 

mes(oj) J u ~ 

The mesh is made out of first order Lagrange triangles, which leads to 2352 degrees 
of freedom. As shown in Figure |T] the structure is decomposed into 9 regular 
subdomains. 
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Figure 1: Substructuring 


We solved the problem with a BDD solver (primal approach). We used the 
Element Equilibration Technique B3 to build statically admissible stress fields; 
each elementary problem was solved using an interpolation of degree 4. For the 
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sake of simplicity, we note: e = ecR n (u N ,Q_ N ) and e = ecR n {u N ,g_N). The results 
in TableQ]are given when the solver is converged (when the residual of both forward 
and adjoint problem are negligible with respect to the loads); they show that the 
error estimation in a substructered context is as efficient as in a sequential context. 
We recall that the slight difference is due to the fact that the static admissible fields 
are not identical in sequential and in domain decomposition methods. 


case 

e 

e 

Ih 

IhH, 2 

iee 

lex 

Sequential 

2.417 

10.23 

2.916 

3.305 10" 4 

12.36 

2.916 

Substructured 

2.418 

10.23 

2.916 

3.858 10~ 4 

12.37 

2.916 


Table 1: Substructured and sequentiel goal-oriented error estimations 


Forward problem Adjoint problem 




Figure 2: Error estimates with separation of sources 


The plots on Figure [2] show the evolution of the two contributions to the bound 
during the iterations of the conjugate gradient. We observe that the discretization 
error is almost constant whereas the algebraic error decreases along the iterations. 
On that example after 6 iterations, the solver errors of both forward and adjoint 
problems are negligible with respect to the contribution of the discretization, extra 
iterations are then useless, they do not improve the quality of the approximation. 
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Improvement of the bound by local refinement 


3.1 Motivation 

The upper bounds obtained in Property [Hand Property [2] always involve the prod¬ 
uct between the error estimators of the forward and adjoint problem. Accordingly 
increasing the quality of either (or both) the forward or the adjoint problems 
necessarily improves the bounds. For instance, since the solution of the adjoint 
problem is highly localized, it is often easier to better its estimation, either by 
mesh refinement or by the injection of a priori knowledge PH HU El E]. 

Adapting a mesh, even for a local quantity of interest, is a global problem: 
simple algorithms like the one based on bisection [3H] propagate on the whole 
domain; parallel remeshing involves exchanges between neighbors (and iterations if 
good load balance is also wished) |£j. Enrichment also implies integrating functions 
on parts of the domain. In this section, we try to benefit the domain decomposition 
in order to simplify these procedure. The idea is to improve the approximation 
only on the subdomains which most contribute to the estimators. We show that a 
class of adaptation algorithms can be realized completely locally without impairing 
the exactness and the computability of the bound: algorithms which modify the 
boundary of the treated subdomains only by adding new degrees of freedom. This 
leads to what we call “nested interfaces” which means that the trace of the initial 
(coarse) kinematic is a subspace of the trace of the new (fine) kinematic. 

Typically, bisections stopped at the boundary of the subdomains belongs to 
that class of algorithms, so do hierarchical h and p refinement and most enrichment 
methods. In the examples, for simplicity reasons, we use uniform h-refinement in 
the subdomains. Of course such procedures open the question of balancing the 
computational load between subdomains, which will be the subject of future work. 
Anyhow, a classical answer is to split the structure into many subdomains and to 
randomly attribute several subdomains to each computing unit with the hope that 
the to-be-refined subdomains will be evenly distributed. 

The rest of the section is organized as follows. In Subsection 13.21 we show how 
a master/slave algorithm enables us to treat nested interfaces in order to preserve 
the ability to build kinematic and static admissible fields. This property would 
be much more difficult to conserve with mortar methods SSi- The algorithm is 
presented in Subsection 13.31 we choose to solve simultaneously the forward and 
the adjoint problems on the locally refined mesh, since there is no real-extra cost 
in solving both problems at the same time. In Subsection 13.41 we verify that the 
chosen strategy enables us to improve the bound on the estimation. 
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3.2 Interpolation of displacements and forces 


We assume that after a first resolution, some subdomains have been locally refined 
(either in h or p sense) or enriched, so that new degrees of freedom have been 
inserted on the boundary of the domain whereas former boundary nodes have 
been preserved. Figure [3] presents a H/2 refinement in the case of 2 subdomains. 
In order to keep the ability of reconstructing globally admissible fields, we decide 
the (coarse) interface to be the master when imposing displacements or forces on 
subdomains, and we keep the assembling operators unchanged. We thus only need 
to derive the interpolation matrix which defines local (fine) displacements 
from interface (coarse) displacement = A^U: = N^u^; and the 

interpolation matrix which defines local (fine) nodal forces A}/ from interface 
) forces A^ = BW t A: A P = M^A^. 


coarse 


1 

1 

1 

Q x 




2 



2 


2 


Interface 


Figure 3: Tit re 


Let <p>P and ipp be the matrices of the coarse and fine shape functions, so 
that displacements write up) = p>p)up' > on the interface and up = cp p uP on the 

boundary. Let be the coordinates of the fine nodes, we have <pP(X.P) = I. 
Imposing the continuity of the displacement writes: 


J s ) = „,«„(*) = 02 (s) U (s) 

—C T_c U C' —F 

^<Pp(X F )^uP = uP 



which defines the interpolation matrix for displacement: 

N (s) = cpp\X.P). (32) 

Matrix N*- 5 ' simply imposes extra fine nodes to follow the displacement imposed 
by the coarse discretization. 

The generation of statically admissible fields requires an algorithm, noted Q 
in [2E], which defines a continuous distribution of traction t from nodal reactions: 
t = Q A. Then transfer operator must ensure that the same continuous 

tractions are created from the coarse and the fine information: t_P = t)p. In other 
words = Q c . 
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Algorithm Q in [28j consists in finding a distribution of traction on the same 
finite element subspace as the displacement, under the constraint to preserve the 
mechanical work: 


ic = ( G c) A g } with G< c 





Xp) with 




(33) 


Gand GP are the L 2 (T) mass matrices of the coarse and fine interpolations. 
Because of the chosen interpolation, it is sufficient to satisfy the relationship = 
t^) on the fine nodes, which leads to: 

M (s) = Gj } N (,) (G^) -1 (34) 


Remark. We also have the following classical expressions for matrices and 


M (s) 


« 


G 


Gfc I G 


G 


< 0 ) 

'c 


(») 

FC 


-1 


► with G^ = 




(35) 


When testing the balance of tractions in the interface displacements we obtain the 
classical preservation of the mechanical work which writes: 


N (s)T M ( s) = I (36) 

Remark. Of course, for unchanged (non-refined) subdomains, we have = 

M<» = I. 


3.3 FETI and BDD algorithms with nonconforming inter¬ 
faces 

The standard BDD algorithm (see for instance Algorithm 1 in [34]) is barely mod¬ 
ified, only the assembling operators change depending if they manipulate displace¬ 
ment or forces: 

• BDD assembling operator becomes A^N^ t , 

• BDD scaled assembling operator becomes A^M^ T , 

The final algorithm is summarized in Algorithm [2] 

The preservation of work (1361) together with the definition of the scaling op¬ 
erators (USD ensure that the reconstructed fields satisfy the required continuity or 
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Algorithm 2: Block-BDD with incompatibility, unknown U = [u, u] 

U = Initialize(F (s) ) ; 

(A^U^) = Solve D (NWAW T U,F( s )); 

Compute residual R = A^N^A^; 

Define local traction SA {s) = M (s )AW t R ; // A^ } = A^ - SA {s) 

<5u (s) = Solves (£A (s >, 0) ; // U$ = Ug } - hU (s) 


Preconditioned residual Z = A^M^^hU^ ; 
Search direction W = PxZ; 
while VR r Z > e do 

(dA^hU^) = Solve D (NWAW T W,0); 

Q = X! s AWnW t <5A£ ) ; 
a = (Q r W)- 1 (R T Z); 

U^U + Wa; 


u£? «- ug } + 

A£? <- A? + 6A$a 


R <— R — Qa; 

(5A (s) = M^A^R ; // A^ } = A^ - (5A (s) 

hU (s) = Solves(<5A (s \ 0) ; // U$ = u£ } - 5U (s) 

Z = ^ s A^M^SXJ^; 

W «- PiZ W(Q T W)- 1 (Q r P 1 Z) 

end 


balance on the interface. More precisely in the BDD algorithm, the displacement 
is built as a continuous field, whereas the nodal forces A$ are balanced on the 
interface because: 


^ A h) N h) T A C 


2 A (s) N (s)T (A^ - 6A {s) ) 

s 

(I — ^ A (s) N (s)T M ( s) A (s)T )R 

s 

(R^A (s) A ( s)T )R = 0 

s 


(37) 


Remark. Similarly, FETI algorithm can be derived from Algorithm 2 in [Mj by 
simply replacing the assembling operators: 

• FETI assembly operator becomes B^M^ T , 

• FETI scaled assembling operator becomes B^N^ T . 
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Remark. When the solver is converged, because of the incompatible meshes, the 
displacement fields (u'd )s and (u$) s are no longer identical. 

3.4 Assessment 

Let us reuse the example of previous section. Figure (Hal) plots the error map of 
the adjoint problem on a first coarse mesh, it is essentially localized within the 
central subdomain; where the loading of the adjoint problem is defined. 


> 0,02 




i 

0 

(a) Coarse mesh (b) Locally refined mesh (c) 

Figure 4: Effect of local refinement (map of the error estimator) on the adjoint 
problem 

Accordingly, from the point of view of the adjoint (resp. forward) problem, we 
need to refine the mesh of the central subdomain. Three cases are considered: two 
with conforming meshes (i) a coarse mesh (given H ) (ii) a fine mesh {h = H/ 4), 
and one with non-conforming meshes due local refinement only in the central sub- 
domain (h = H/4, see Figure [5] for a schematic presentation with a coarser mesh 
than in actual computations). In all cases, we compute the error estimator at con¬ 
vergence of the solver. The number of Conjugate Gradient iterations is identical 
for the 3 meshes (9 iterations). The results are summarized in Table [21 When the 
central subdomain is refined, the error estimator for the adjoint problem decreases 
while the error estimator of the forward problem remains more or less unchanged. 
As a consequence, the use of local refinement (independently by subdomains) en¬ 
ables us to obtain improved bounds on the quantity of interest (—32%) with a 
simple implementation and less degrees of freedom (+38% #dof) than would be 
required by the full fine mesh (+300% ^dof for a —65% decrease of the estimator). 

Figure |4b] plots the map of error for the adjoint problem when only the central 
subdomain is refined. The source of error due to the non-conforming interface 
(boundary of the central subdomain) is clearly negligible. 


18 





Figure 5: Hierarchically refined mesh on the central subdomain 


case 

e 

e 

Ih 

IhH, 2 

¥*■ 

#dof 

Coarse mesh 
Fine mesh 

2.418 

1.216 

10.23 

7.020 

2.916 

2.916 

3.8575 10~ 4 
-9.695 10~ 4 

12.37 

4.268 

8738 

34988 

Refined central SD 

2.404 

7.020 

2.9169 

3.8575 10” 4 

8.438 

12122 


Table 2: Comparaison between conforming and non conforming cases, I ex = 2, 916 


Of course, the use of locally refined meshes is questionable if the zone of high 
contribution to the error is close to a non-matching interface: in that case the 
driving of the fine kinematic by the coarse interface is an extra source of error. 
To evaluate that problem, we study how the efficiency of the local refinement is 
altered by the distance d between the support of the quantity of interest u = 
[— yq — d] jq — d] x [— {q] j^] and the interface T. 

For each configuration (characterized by the distance d) we compute the relative 
discretization error (discretization error divided by the energy of the finite element 
solution) of the adjoint problem p^h (resp. ph) when h = H/ 4 only in the central 
subdomain (resp. on all subdomains). Figure [6] plots the ratio pHh/ Ph as a function 
of the distance d/do (where d = do corresponds to the quantity of interest at the 
center of the subdomain). The quality of the refined estimator is degraded only 
when d/do < 10” 1 which corresponds to a quantity of interest very close to the 
interface. Thus the idea to conduct local refinements remains valid in a wide range 
of configurations. 
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do 


Figure 6: Influence of the distance d 


4 Numerical experiments 


The aim of these assessments is to show how easily the error bound can be improved 
by the refinement strategy described in Section [3j 

We always use the stopping criterion proposed in Section [2j iterations are 
stopped when the algebraic error becomes negligible (10 times smaller) with respect 
to the contribution of the discretization estimated at the first iteration, for both 
forward and adjoint problems. We verify that the discretization error computed 
with the new interface fields barely varies, as could be expected from [53] . 

After a first resolution on a coarse mesh, we analyze the contributions of sub- 
domains to the forward and adjoint error estimator using the following measures: 


= 


e CR n(s) 


'CR, 


n( s ) 


\U 


(d £(») 


and 


^0) = 




0) 
■N i 
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(W 

■N / 


■N > 


■N 


fn (s) 


^ e? 


(38) 


From that information we deduce h-refinement strategies and evaluate their ability 
to improve the estimator. 


4.1 The Gamma-shaped structure 

We consider a Gamma-shaped structure clamped on its base and submitted to 
prescribed traction and shear on the upper-right side. The problem is linear elastic 
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with unit Young modulus and 0.3 Poisson coefficient. The initial mesh is made out 
of first order triangles. The structure is decomposed into 7 domains and we use 
the classical BDD solver. The quantity of interest is the mean value of the strain 
e yy over the region oj located in the seventh subdomain, as illustrated in Figure 
[7] We build statically admissible stress fields with the Flux-free technique [71 [27]. 
Each star-patch problem is solved with a forth order interpolation. 



Figure 7: Loading of forward (blue) and adjoint problems (orange) 

The contributions 7p s ' and rf 8 ^ to the error estimator on the initial mesh are 
presented in Figure [8aJ We observe that the major contribution to the error for 
the forward problem is located in the fifth subdomain (because of the corner). The 
objective is to refine this domain in order to obtain a level of error comparable 
to subdomain 1 (which is the second largest contributor). It implies to divide the 
error in the hfth subdomain by 2.8. The convergence in this subdomain is driven 
by the corner singularity. Therefore, we decide to perform three hierarchical h- 
reffiiements in the hfth subdomain. The quantity of interest being localized in 
Subdomain 7, we also process 3 hierarchical h-rehnements in that subdomain. 

As expected, we can see on Figure [8b] that the forward error in Subdomain 5 
and Subdomain 1 are nearly the same on the new mesh. We also observe that the 
error of the adjoint problem has decreased (due to the refinement of Subdomain 7). 
In Table [3] we sum up the global values of discretization errors, the quantities Ih 
and Ihh, 2 - 


Mesh 

e 

e 

Ih 

IhH, 2 

iee 

Uniform 

33.528 

0.20977 

-0.39157 

-0.0010045 

3.5165 

Locally refined 

18.464 

0.072916 

-0.39006 

-0.00050892 

0.67316 


Table 3: Performance of local refinement on the T-shaped structure 
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(a) Initial mesh 
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0.8 - e = 18.464 
e = 0.072916 

0.6 - 


0.4 - 

0.2 - 

0 - 


1 2 3 4 5 6 

(b) Locally refined mesh 


Figure 8: Distribution of the error estimator within subdomains 


4.2 The pre-cracked structure 

We now consider the pre-cracked structure decomposed into 16 subdomains as 
illustrated in Figure [9j The displacements at the base of the structure and on the 
bigger hole are imposed to be zero. The upper-left part and the second hole are 
subjected to a constant unit pressure. The quantity of interest is the mean of the 
stress component a xx on a region oj close to the crack. In Figure |H1 the loading of 
the reference problem is in blue and the loading of adjoint problem is in orange. 
We used the FETI algorithm (dual approach) to solve the interface problem and 
the statically admissible stress fields are built using the Flux-free technique Hi El. 
with forth order interpolation for problems on star-patches. 
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(a) Initial mesh (b) Locally refined mesh 


Figure 10: Distribution of error within subdomains 


In Figure IIPal we observe that the discretization errors of the forward and 
adjoint problem are quasi totally located in the sixth subdomain. Therefore, we 
perform hierarchical refinement in this subdomain in order to reach a level of 
forward error comparable to the second larger contributor (Subdomain 8). That 
corresponds to dividing the error in Subdomain 6 by a factor 3. Because of the 
singularity in the subdomain, convergence is expected to evolve like \[h so we 
perform 4 hierarchical refinements (each coarse element edge is divided into 16 
fine edges). 

As can be seen in Figure [TObl after local refinement the contributions of Sub- 
domains 6 and 8 to the forward error are equivalent, and the adjoint error also 
decreased, as expected. The global performance is summed up in Table [4] In the 
end, the bound on the estimator on the quantity of interest was decreased by a 
factor 3. 


Mesh 

e 

e 

Ih 

IhH, 2 

iee 

Uniform 

26.215 

0.98905 

2.4915 

3.1935 

12.964 

Locally refined 

16.662 

0.51378 

3.2165 

0.086055 

4.2803 


Table 4: Performance of local refined for the cracked structure 


5 Conclusion 

In this article, we extended two classical results of goal-oriented error estimation 
to problems solved by domain decomposition methods. We solved simultaneously 
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both forward and adjoint problems on the same mesli using a block conjugate 
gradient algorithm and proposed a fully parallel procedure to recover admissible 
fields and estimate the error on the quantity of interest. Moreover, the error 
bounds separated the contribution of the discretization from the contribution of 
the iterative solver which lead to defining of a stopping criterion of the iterative 
solver that prevented oversolving. 

We also showed that the method could be easily extended to the case of locally 
enriched kinematics (using h or p refinements or dedicated enrichments) only con¬ 
strained to contain the original coarse kinematic on the boundary of subdomains. 
This property provided a simple parallel way to improve estimators as soon as 
the singularities (local quantity of interests, crack tips, etc) were not too close to 
the interface. The concepts were validated on academic examples using the most 
basic /(-refinement strategy. Of course using ad-hoc enrichment instead of brutal 
/(-refinement should improve the performance of the method. 

Future work will focus on exploiting the fact that the shape of the interface 
is preserved by the enrichment to define acceleration strategies, and on defining 
original strategies to balance the computational load between subdomains for the 
resolution of the refined problem. 
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